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ABSTRACT 


The Eurasian ice sheet complex (EISC) was the third largest ice mass during the Last Glacial Maximum 
with a span of over 4500 km and responsible for around 20 m of eustatic sea-level lowering. Whilst 
recent terrestrial and marine empirical insights have improved understanding of the chronology, pattern 
and rates of retreat of this vast ice sheet, a concerted attempt to model the deglaciation of the EISC 
honouring these new constraints is conspicuously lacking. Here, we apply a first-order, thermo- 
mechanical ice sheet model, validated against a diverse suite of empirical data, to investigate the retreat 
of the EISC after 23 ka BP, directly extending the work of Patton et al. (2016) who modelled the build-up 
to its maximum extent. Retreat of the ice sheet complex was highly asynchronous, reflecting contrasting 
regional sensitivities to climate forcing, oceanic influence, and internal dynamics. Most rapid retreat was 
experienced across the Barents Sea sector after 17.8 ka BP when this marine-based ice sheet disintegrated 
at a rate of ~670 gigatonnes per year (Gt a`!) through enhanced calving and interior dynamic thinning, 
driven by oceanic/atmospheric warming and exacerbated by eustatic sea-level rise. From 14.9 to 12.9 ka 
BP the EISC lost on average 750 Gt a`!, peaking at rates >3000 Gt a~', roughly equally partitioned be- 
tween surface melt and dynamic losses, and potentially contributing up to 2.5 m to global sea-level rise 
during Meltwater Pulse 1A. Independent glacio-isostatic modelling constrained by an extensive in- 
ventory of relative sea-level change corroborates our ice sheet loading history of the Barents Sea sector. 
Subglacial conditions were predominately temperate during deglaciation, with over 6000 subglacial 
lakes predicted along with an extensive subglacial drainage network. Moreover, the maximum EISC and 
its isostatic footprint had a profound impact on the proglacial hydrological network, forming the Fleuve 
Manche mega-catchment which had an area of ~2.5 x 10° km? and drained the present day Vistula, Elbe, 
Rhine and Thames rivers through the Seine Estuary. During the Bglling/Allerad oscillation after c. 14.6 ka 
BP, two major proglacial lakes formed in the Baltic and White seas, buffering meltwater pulses from 
eastern Fennoscandia through to the Younger Dryas when these massive proglacial freshwater lakes 
flooded into the North Atlantic Ocean. Deglaciation temporarily abated during the Younger Dryas stadial 
at 12.9 ka BP, when remnant ice across Svalbard, Franz Josef Land, Novaya Zemlya, Fennoscandia and 
Scotland experienced a short-lived but dynamic re-advance. The final stage of deglaciation converged on 
present day ice cover around the Scandes mountains and the Barents Sea by 8.7 ka BP, although the 

phase-lagged isostatic recovery still continues today. 
© 2017 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY-NC-ND 
license (http://creativecommons.org/licenses/by-nc-nd/4.0/). 
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1. Introduction 


Northern Eurasia was covered by three semi-independent ice 
sheets that between 26 and 19 ka BP (Clark et al., 2009) coalesced to 
form a single Eurasian ice sheet complex (EISC) during the last 
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glacial maximum (LGM) (Svendsen et al., 2004). This complex had 
an impressive latitudinal and longitudinal coverage, with contin- 
uous ice cover spanning over 4500 km extending southwest of the 
Isles of Scilly (50°N, 6°W) on the Atlantic seaboard to beyond Franz 
Josef Land (81°N, 56°E) in the Russian High Arctic (Fig. 1). It was the 
third largest ice mass after the North American and Antarctic ice 
sheets, and, with a combined volume three times the present 
Greenland ice sheet, accounted for at least 20 m of eustatic sea- 
level lowering (Patton et al., 2016). Growth of the EISC initiated 
from three main nucleation centres located over Britain and 
Ireland, Fennoscandia, and the Barents-Kara seas, with contrasting 
styles of glaciation and associated conditions and processes 
reflecting these settings, from marine-terminating, fast-flowing ice 
streams in maritime regions to extensive frozen-based glaciation in 
inter-ice-stream and upland areas. 

Knowledge of the maximum extent, chronology, and patterns of 
retreat of the EISC has improved greatly in the last decade, partic- 
ularly in offshore sectors where marine geophysical surveys have 
addressed a number of notable gaps in understanding (e.g., Landvik 
et al., 2005; Ottesen et al., 2005; Bradwell et al., 2008; Dunlop et al., 
2010; Winsborrow et al., 2010; Andreassen et al., 2014; Sejrup et al., 
2016). Moreover, developments in cosmogenic exposure dating, as 
well as refinement of radiocarbon and other dating techniques, has 
enabled detailed onshore deglaciation chronologies to be devel- 
oped (Rinterknecht et al., 2006; Linge et al., 2007; Ballantyne, 2010; 
Stroeven et al., 2011; Briner et al., 2016). Subsequent publication of 
data-rich compilations and review studies has thus set in place a 
strengthened empirical framework against which modelling in- 
vestigations of the EISC can be made and tested (Napieralski et al., 
2007; Clark et al., 2012; Hormes et al., 2013; Hughes et al., 2014, 
2016; Patton et al., 2015; Auriac et al., 2016; Cuzzone et al., 2016; 
Stroeven et al., 2016), although the utility of such empirical data- 
sets for modelled reconstruction comparisons must be considered 
in light of the quality of legacy data they incorporate (e.g., Small 
et al., 2017). 

Since the early numerical modelling undertaken as part of the 
QUEEN programme (cf. Siegert and Dowdeswell, 2004), progress on 
modelling the Late Glacial retreat of the EISC has been limited to a 
number of regional reconstructions (Holmlund and Fastook, 1993; 
Boulton et al., 2001; Boulton and Hagdorn, 2006; Hubbard et al., 
2009) or otherwise focussed primarily on process dynamics 
(Arnold and Sharp, 2002; Forsstrom and Greve, 2004; Naslund 
et al., 2005; van den Berg et al., 2008a, 2008b, Clason et al., 2014, 
2016). An alternative to these process-based models are EISC re- 
constructions developed through glacial isostatic adjustment 
modelling (Peltier, 2004; Lambeck et al., 2010; Peltier et al., 2015). 
These inverse models are calibrated using empirically determined 
ice extents and relative sea-level data, but the resulting re- 
constructions are static and do not provide insight into the dy- 
namics of ice sheet retreat, nor do they inform the climatic/oceanic 
forcing that drove it. 

In this paper, we apply a first-order, thermomechanical ice sheet 
model to investigate the dynamic retreat of the EISC after 23 ka BP. 
The primary aims are twofold: i) to present a robust, 4D high- 
resolution, synoptic reconstruction of EISC deglaciation from 23 
to 8 ka BP, from its local LGM extent, through the Younger Dryas 
stadial (12.9—11.7 ka BP), and into the Early Holocene; and, ii) to 
validate and discuss model output against a suite of empirical data 
that constrain both the pattern and rate of retreat of the EISC, 
including its glacial-isostatic footprint, chronological data for the 
timing of deglaciation, flowset vectors, and its sub- and pro-glacial 
hydrological legacy. 

The study extends the work of Patton et al. (2016) who previ- 
ously explored the asynchronous and asymmetric growth of the 
EISC to its maximum LGM extent from 37 to 19 ka BP. A variety of 


geomorphological, geophysical and geochronological data are used 
to constrain and validate the broad-scale dynamics of the retreating 
ice mass. In particular though, where terrestrial constraints are 
notably lacking across the relatively data-sparse Barents and Kara 
seas, we utilise glacial isostatic adjustment modelling to test tran- 
sient ice loading and retreat rates. 


2. Methods 
2.1. The ice flow model 


The 3D thermomechanical model and associated initial bound- 
ary condition data applied are of the same derivation as that used to 
previously model the pre-LGM build-up of the EISC by Patton et al. 
(2016), where a more complete description of the model setup and 
implementation can be found. In brief, the ice flow model is a first- 
order approximation of the Stokes equations, adapted from Blatter 
(1995), Hubbard (1999, 2000), Marshall et al. (2005), and Pollard 
and DeConto (2007). The approach to solving the three dimen- 
sional stress/strain field equates to the L1L2 classification of higher- 
order models defined by Hindmarsh (2004), and includes longitu- 
dinal (membrane) stresses that become increasingly important 
across steep gradients in topography and motion. The model is 
integrated forward through time on a finite-difference grid with a 
resolution of 10 km through perturbations in climate and eustatic 
sea level (Fig. 2A—B). Isostatic loading is implemented using an 
elastic lithosphere/relaxed asthenosphere scheme described by Le 
Meur and Huybrechts (1996), which provides a computationally 
pragmatic solution in the absence of a full spherical earth model. 
Gridded output is projected under an equal area Lambert 
Azimuthal projection, with a central meridian of 73°E. 

Surface mass balance is determined by a positive degree-day 
scheme, applied according to Laumann and Reeh (1993), and de- 
rives total melt from integrated monthly positive temperatures. 
Both temperature and precipitation adjust to the evolving ice sheet 
surface through applied lapse rates derived from multiple- 
regression analyses of meteorological observations at a resolution 
of 1 km from the WorldClim database (Hijmans et al., 2005; Version 
1.4). To account for the large variations in climate regime across the 
Eurasian domain, regional reference climates and associated forc- 
ing are tuned independently for each of the three major accumu- 
lation centres (Fig. 2C—D). An additional mass balance term 
incorporated is the net water vapour flux to and from the ice sheet 
surface — a predominant component of ablation in cold continental 
settings where humidity can be very low (e.g., Fujii and Kusunoki, 
1982; Kameda et al., 1997). 

Calving losses at marine-terminating margins are coupled to 
relative sea level (RSL) (Waelbroeck et al., 2002) using a standard 
empirical function relating the calving flux to ice thickness and 
water depth (Brown et al., 1982; van der Veen, 1999). The sensitivity 
of calving to, for example, variations in ocean temperature 
(Luckman et al., 2015) and sea-ice buttressing (Hoff et al., 2016) is 
controlled spatially and temporally through a depth-scaled calving 
parameterisation (Hubbard, 2006) (Fig. 2D). In the absence of 
explicit calculations of such external feedbacks, this depth-related 
calving coefficient provides a pragmatic and computationally effi- 
cient parameterisation for determining mass loss at marine ter- 
minating margins of the ice complex. The model is applied to a 
10 km finite-difference mesh with the inclusion of grounding-line 
dynamics based on the analytical boundary-treatment of Schoof 
(2007) and adapted in 2D by Pollard and DeConto (2007), which 
defines the ice flux at the grounding line as a function of ice 
thickness linearly interpolated between the adjacent node that 
bracket floating and grounded ice (Hubbard et al., 2009). 
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Fig. 1. Major drainage routes of the Eurasian ice sheet complex, adapted from Stokes and Clark (2001), Ottesen et al. (2005) and Clark et al. (2012). Locations of major trough mouth 
fans (brown) adapted from Dahlgren et al. (2005) and Batchelor and Dowdeswell (2014). PB: Porcupine Bank; BDF: Barra and Donegal Fans; RB: Rosemary Bank; NSF: North Sea Fan; 
Bj: Bjorngyrenna Fan. Glacial limits are compiled from Ehlers and Gibbard (2007), Patton et al. (2015) and Stroeven et al. (2016). Recent evidence for the extension of the Celtic ice 


sheet onto Porcupine Bank (PB) (Peters et al., 2015) and into the southern Celtic Sea (Praeg et al., 2015) is also incorporated. (For interpretation of the references to colour in this 
figure legend, the reader is referred to the web version of this article.) 
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Fig. 2. Forcing curves used to drive the evolution of the optimum EISC deglaciation experiment (23-8 ka BP), with bulk changes in climate forcing and calving sensitivity tuned 
separately for the three sub-domains of the complex. A) Eustatic sea level forcing (Waelbroeck et al., 2002); B) Mean annual air temperature (MAAT) forcing curves used to drive the 
deglaciation experiments, scaled against the NGRIP 5!80 record (Andersen et al., 2004), C) Transient changes in total precipitation (relative to present day levels) across the three sub 
domains of the ice complex; D) Transient changes in the sensitivity to calving losses, with higher values leading to increased calving losses along a linear relationship. Arctic areas 
are defined as being above 78°N. Chronozones marked include the Bølling/Allerød (Bø-Al) oscillation (14.6—12.9 ka BP) and Younger Dryas (YD) stadial (12.9—11.7 ka BP). 


2.2. Experiment setup and empirical constraints 


The deglaciation scenarios presented in this study follow on 
directly from the suite of EISC build-up experiments carried out by 
Patton et al. (2016), and thus forms part of a continuous recon- 
struction for the full Late Weichselian (37-8 ka BP). With a greater 
abundance of empirical evidence constraining the deglaciation of 
the ice sheet, the prescribing of model parameters in order to yield 
a transient experiment that follows the broad pattern of ice retreat 
becomes increasingly challenging. This challenge is greatly aided by 
the implementation of three sub-domains that span each of the 
major nucleation centres of the complex, within which climatic and 
oceanographic forcings can be tuned independently. During the 
reconstructed build-up of the ice complex, major climatic differ- 
ences between these three regions included a general enhance- 
ment of climate cooling and a decrease in precipitation across 
southern latitudes. To encourage growth of the marine-based 
Barents Sea ice sheet (BSIS), sensitivity to calving losses was also 
reduced relative to that of the Fennoscandian and Celtic ice sheets, 
simulating the influence of colder ocean temperatures and but- 
tressing effects from perennial sea ice (Patton et al., 2016). 

Various lines of geological and geophysical data were used to 
guide the broad retreat patterns of the modelled ice sheet, 
including key dating constraints, geomorphological flow sets, 
notable moraine/grounding zone wedge positions, and patterns of 
isostatic loading. The most important datasets in this respect are 
those that incorporate a holistic approach to the ice-complex 


reconstruction, for example, the DATED-1 1 ka-interval timeslice 
reconstructions made using an extensive database of chronological 
constraints (Hughes et al., 2016). The value of this comprehensive 
approach becomes apparent where conflicting or competing 
regional interpretations occur, as the data can be rigorously 
examined within a broader and glaciologically consistent context, 
and thus provide a sense of probability or even highlight unreliable 
data points. To this end, ice sheet modelling fulfils a useful and 
independent tool for comparison, exploring the feasibility and 
probability of competing hypotheses for ice evolution and dynamic 
behaviour. 

A useful distinction can be made between the pattern and 
timing of deglaciation in constraining and testing the deglaciation 
scenario. Except for the Barents and Kara seas, coverage of 
geomorphological data has expanded hugely in recent years 
through interpretation of terrestrial LIDAR elevation data and sea- 
floor multibeam mapping (e.g., Winsborrow et al., 2010; Hughes 
et al., 2014; Greenwood et al., 2015). This wealth of data provides 
assurances that the broad pattern of deglaciation of the EISC has 
been identified correctly, including the main stillstands and read- 
vances. However, the time sequence of deglaciation of the EISC 
remains less certain. This deficit is partly a product of the limited 
availability of dating controls for Russian territories and the floor of 
the Arctic seas. Yet existing dates elsewhere currently often lack 
precision at millennial time scales. A comprehensive compilation of 
radiocarbon, optically-stimulated luminescence and '°Be cosmo- 
genic isotope ages for the last Fennoscandian ice sheet (FIS) has 
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shown differences of several thousand years in dates for major 
stillstands, underlining the greater precision of the established 
varve chronology (Stroeven et al., 2016). Large Be datasets for 
individual moraines (Rinterknecht et al., 2006) and for specific sites 
(Ballantyne, 2010) also show wide date ranges. Maximising the 
utility of legacy geochronological data for the purposes of inform- 
ing numerical modelling studies has only recently come to the fore 
following the compilation of datasets that span entire ice sheet 
domains (e.g., Dyke et al., 2002; Hughes et al., 2011, 2016; Bentley 
et al., 2014; Stroeven et al., 2016), driving standards for sampling 
strategies and criteria for data reporting (e.g., Small et al., 2017). In 
this respect, the current chronology of EISC deglaciation is insuffi- 
ciently robust, particularly early in deglaciation, to pin down the 
timings of millennial-scale events. Indeed part of the value of this 
simulation lies in the independent constraints that it provides on 
the likely timings of major events during EISC decay. An assessment 
on the conformity of the optimum deglaciation scenario against the 
empirical datasets used is presented in section 4, following the 
description of the modelled deglaciation. 


2.3. Glacial isostatic adjustment modelling 


Ice thickness predictions derived from the ice model above were 
incorporated into an independent glacial isostatic adjustment 
model to examine the overall fit of ice loading to sea-level histories 
across the Barents Sea. The sea-level equation (first derived by 
Farrell and Clark (1976)) is solved using the implementation from 
Mitrovica and Milne (2003) and Kendall et al. (2005). Gravitation- 
ally self-consistent sea-level changes are computed, taking into 
account shoreline evolution as well as the time-dependent evolu- 
tion of marine-based ice margins. The sea-level equation is solved 
iteratively using an extended pseudo-spectral algorithm. This nu- 
merical code assumes a spherically symmetric Earth, whose prop- 
erties are based on the Preliminary Reference Earth Model (PREM; 
Dziewonski and Anderson, 1981). The Earth model is implemented 
as an input with three variables: lithosphere thickness and upper 
and lower mantle viscosity. We use 300 different Earth models, 
where the lithosphere thickness ranges from 46 to 120 km and the 
upper and lower mantle viscosities range from 0.05 x 102! to 
5 x 10°! Pasand 1 x 107! to 50 x 107! Pas, respectively. These Earth 
models cover the range of Earth parameters generally found or 
inferred for this area from a range of geophysical techniques (e.g., 
Kaufmann and Wolf, 1996; Steffen and Kaufmann, 2005; Klitzke 
et al., 2015). After solving the sea-level equation, an estimate of 
the time evolution of relative sea level at specific locations is 
determined. 


3. Results 
3.1. Ice sheet build-up 


The deglaciation scenario presented in this paper forms a direct 
continuation of the optimum build-up experiment (37-19 ka BP) 
presented by Patton et al. (2016). The progression and dynamic 
behaviour of the retreating ice complex is therefore, at least 
initially, predetermined by this pre-LGM growth trajectory. For 
example, the eastwards migration of the central ice divide of the 
complex that forces asynchronous maximum extension between 
25 and 20 ka BP plays a dominant role in the subsequent timing and 
direction of ice retreat. In order to place the new experiment results 
in context, a brief description of the modelled build-up of the EISC 
by is provided here; areal and volumetric changes for each semi- 
independent ice sheet through the full Late Weichselian are pro- 
vided in Fig. 3. 

From a distribution of limited ice cover, widespread ice growth 


initiates across all three sectors in response to climate deterioration 
during the latter stages of marine oxygen isotope stage (MIS) 3 (<37 
ka BP). Continuous net accumulation sees the emergence of fast 
flowing outlet glaciers that reach the shelf break west of Svalbard, 
Norway, Britain, and Ireland around 30-29 ka BP, indicative of 
extensive and thick terrestrial ice sheets. By 24.5 ka BP the 
expanding Barents Sea and Fennoscandia ice sheets coalesce, 
forming a major ice divide running northwards over Finnmark and 
the central banks of the Barents Sea. The Celtic and Fennoscandia 
ice sheets merge in the North Sea 1.5 ka later, completing the for- 
mation of a coherent Eurasian ice sheet complex by 23 ka BP. With 
stable western margins between 25 and 23 ka BP, continuing ice 
divide migrations force an ice incursion into eastern European 
sectors as late as 20 ka BP. A weak coalescence between the Celtic 
and Fennoscandian ice sheets sees the two ice sheets split shortly 
after 22.4 ka BP, with the major Norwegian Channel ice stream 
persisting in the North Sea until at least 19 ka BP. Shortcomings of 
this build-up scenario that impinge on the accuracy of the degla- 
ciation scenario presented include the limited expansion of the 
Norwegian Channel ice stream ~400 km short from the continental 
shelf edge at 22.7 ka BP, the failure to reproduce the extended 
outlet lobes of the eastern FIS (Dvina, Vologda and Rybinsk basins), 
and the shortfall of the maximum modelled Fennoscandian ice 
margin in Germany and Poland. 


3.2. Eurasian ice sheet complex deglaciation 


3.2.1. 23.0—17.8 ka BP 

A period of climate stability c. 23.0—19.5 ka BP (Fig. 2B) sees 
relatively minor volume and areal fluctuations across the entire 
complex (Fig. 3). The furthest extent of the EISC through the White 
Sea and into northwest Russia occurs early in the deglaciation 
sequence at 20.4 ka BP (Fig. 4A). The asynchronous development of 
the LGM complex resonates with empirical observations (Hughes 
et al., 2016; Stroeven et al., 2016), with maximum extension of 
the north-eastern sectors of the FIS possibly occurring as late as 18- 
16 ka BP (Larsen et al., 1999b, 2016). Elsewhere, retreat of westward 
draining ice along much of the western continental shelf edge is 
interrupted, stabilising at positions upslope from the shelf break. A 
phase of warming after 19.5 ka BP is most severely felt across Celtic 
sectors, the North Sea Basin, and the south eastern FIS. Ice retreat 
occurs throughout much of Ireland, Wales, northern England, with 
substantial ice drainage occurring through the Malin Sea, and 
southern Scandinavia. However, the pace of retreat is stepwise, 
with a brief respite linked to mean annual air temperature (MAAT) 
cooling at 17.8 ka BP (Fig. 4B). Surface melting across Fennoscandia 
and the Barents Sea at this time has a limited impact on ice extent, 
with major outlet glaciers undergoing minimal retreat. 

A manually forced reduction of the calving coefficient along the 
Arctic Ocean margin (>78°N), resulting in negligible calving losses, 
enforces a largely stationary margin throughout this period, and 
was found to be a necessary forcing in order to satisfy relative sea- 
level constraints on Franz Josef Land and Novaya Zemlya (cf. section 
3.3). However, fluctuations in MAAT force at least three phases of 
stillstand/readvance within Bjgrngyrenna (Fig. 4B) before 17.8 ka 
BP, some 100 km away from the shelf break. 


3.2.2. 17.8—14.9 ka BP 

The minor standstill at 17.8 ka BP for some margins of EISC is 
interrupted by a phase of rapid retreat at 17.5 ka BP, for example, 
through Bjgrngyrenna and the coast parallel troughs of northern 
Fennoscandia (Fig. 4C), with only a minor re-advance c. 16 ka BP 
interrupting continuous retreat (Fig. 4D). Vigorous draw-down of 
ice through the various troughs of the central Barents Sea forces an 
effective collapse of this ice saddle, with the BSIS and FIS eventually 
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Fig. 3. Areal and volumetric changes for the three semi-independent ice sheet components of the EISC through the full Late Weichselian, extending the optimum build-up (37—19 
ka BP) experiment of Patton et al. (2016). Acronyms used: CIS (Celtic ice sheet); FIS (Fennoscandian ice sheet); and BSIS (Barents Sea ice sheet). Note the log scale for CIS ice volume. 
Chronozones marked include the Bglling/Allerad (Bg-Al) oscillation (14.6—12.9 ka BP) and Younger Dryas (YD) stadial (12.9—11.7 ka BP). 


separating at 15.5 ka BP (Fig. 4D—E). Chronological control on this 
separation is poor, and is constrained here by a single radiocarbon 
date of 14.9 cal ka BP from glaciomarine sediments in the outer 
Pechora Basin (Polyak et al., 1995). Further calving losses continue 
to force rapid withdrawal of the ice margin northwards through the 
Barents Sea, with little topographic relief available to provide stable 
pinning points. 

Despite widespread ice losses from all margins of the FIS after 
17.5 ka BP, the ice sheet stabilises c. 1 ka later, with various outlet 
glaciers, including those in the Norwegian Channel, Baltic Sea and 
White Sea, even undergoing a phase of readvance after c. 16 ka BP 
until the end of Greenland Stadial-2.1a at 14.7 ka BP (Fig. 4E). 

For the diminishing, yet sensitive, Celtic ice sheet (CIS), the Late 
Glacial warming is critical. By 16.5 ka BP the CIS is all but pinned to 
the Scottish coastline until the end of GS-2.1a. During this period, 
active drainage routes include the Moray Firth, Minch, and Sea of 
the Hebrides (Fig. 4E). 


3.2.3. 14.9—12.9 ka BP 

The climate warming that marks the start of the Belling-Allerad 
oscillation inflicts a period of widespread melt, at times with mass 
losses >3000 Gt a~! occurring across the entire complex (Fig. 3B). 
During this period (14.9—12.9 ka BP) the EISC loses on average 
750 Gt a, partitioned roughly equally between surface melt and 
dynamic losses, with the EISC contribution to Meltwater Pulse 1A 
(Deschamps et al., 2012) on the order of ~2.5 m of sea-level 
equivalent ice volume. Within 500 years of the initiation of the 
Bølling, the remnants of the CIS are but completely disappeared 
(Fig. 4F). Combined with continued calving losses, this period is 
equally catastrophic for the BSIS — by the start of the Younger Dryas 
at 12.9 ka BP all that is left of the ice sheet are isolated ice caps 
across the islands, and an ice saddle between eastern Svalbard and 
Franz Josef Land (Fig. 4G). For the FIS the greatest losses occur along 
the southern and eastern margins furthest away from the central 
ice divide, forcing a westward migration of the ice divide. Although 
the western ice margin continues to be pinned along the Norwe- 
gian coast throughout this period, retreat of ice through Sweden, 
Finland and northwest Russia continues unabated until 12.9 ka BP, 


when substantial drops in MAAT halt retreat. 

Although proto versions of the Baltic and White Sea proglacial 
lakes are in place prior to the Balling-Allerad oscillation (Fig. 4B—E), 
their capacity only dramatically expands with the onset of ice 
retreat through the respective basins during this interstadial, as 
suggested by Björck (2008). 


3.2.4. 12.9—8.0 ka BP 

The initiation of the Younger Dryas at 12.9 ka BP sees a re- 
emergence of cold-based ice cover over the Scottish uplands, 
though this is short-lived and disappears soon after the end of the 
Younger Dryas by 11.4 ka BP (Fig. 4G). The response of the FIS is 
somewhat different, with the Bothnian Sea ice stream continuing to 
retreat through the beginning of this stadial despite thickening of 
ice at the divide. The FIS eventually purges mass at 12.3 ka BP, 
undergoing a rapid readvance mainly across its southern margin, 
but also offshore into the Norwegian Sea. The readvance culminates 
at 11.7 ka BP, after which warming during the early Holocene forces 
terminal retreat of ice into two separate ice caps across the Scandes. 
The smaller of the two is focussed over the southern mountains 
around Jotunheimen, while the larger retreats to the northern 
Swedish mountains (Fig. 4H). The death of the Bothnian Sea ice 
stream is marked by a phase of rapid flow with an onset zone close 
to the present-day coastline at c. 10 ka BP. The northern FIS centre 
finally disappears at 8.7 ka BP. 

Despite the climate deterioration during the Younger Dryas, 
marine-based sectors of the BSIS finally disappear during this 
period, with the ice saddle connecting Svalbard and Franz Josef 
Land collapsing at 12.2 ka BP. Reduced moisture supply at this time 
is insufficient to drive a readvance of the remaining ice caps, and 
warming through the Holocene continues the overall trend of 
margin retreat. 


3.3. Glacial isostatic adjustment modelling 
Results from the comparison between our optimum deglacial 


scenario and 46 (RSL) observations from around the Barents Sea are 
given in Table 1. These are subdivided into the four main terrestrial 
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Fig. 4. Timeslices showing ice extent and ice-surface velocities during deglaciation of the EISC between 20.4—10.0 ka BP. Surface elevation contours are drawn every 250 m. Sea 
level is determined from the global sea level forcing curve (Fig. 2A) and includes corrections based on modelled isostatic loading effects. Shorelines of major proglacial lakes are 
defined based on the elevation of the contemporary sea level position. Acronyms used: G — Goteborg Moraine; HC — Halland Coastal Moraine; SoH — Sea of Hebrides. 


areas bordering the domain: Svalbard, Franz Josef Land, Novaya 
Zemlya and northern Fennoscandia. RSL curves, selected as being 
representative of the full array of results are presented in Fig. 5. The 
full set of RSL plots can be found as supplementary material in 
Figure S1. Table 1 presents the best-fit earth model parameters for 
each region of the Barents Sea, as well as the corresponding value of 
x2 — a measure of data fit derived from the mean weighted residual 
sum of squares between the modelled and empirical data points. 
These values are compared to an earlier published iteration of the 


modelled deglaciation wherein ice retreat was more loosely con- 
strained (Auriac et al., 2016). 

While the global best fit Oe ) is still well above 1, indicating that 
the ice load scenario fails to reproduce the RSL observations 
simultaneously at all sites, work to better constrain the 3D 
geometrical evolution of the BSIS has yielded improvements if 
regional variations to the Earth model are assumed (Table 1). 
Northern Fennoscandia shows the most clear improvement, with a 
decrease in overall misfit (x2) from 51.9 to 8.13. Across Franz Josef 


H. Patton et al. / Quaternary Science Reviews 169 (2017) 148—172 155 


Table 1 


x2 and best-fitting Earth rheology parameters for regional fits of the optimum reconstruction, compared with an earlier and less well-constrained model experiment 
(UiT_2016) of the deglaciation presented by Auriac et al. (2016). Parameters include: h (lithosphere thickness — km), v, (upper mantle viscosity—10! Pa s), and v; (lower mantle 


viscosity—107! Pa s). 


Model Svalbard (n = 18) Franz Josef Land (n = 13) 
XZ h Vu v 2 h vu 
UiT_2016 20.6 120 0.2 2 2.5 96 3 
This paper 75.3 120 1 1 31 120 3 
Model Svalbard (East)* (n = 11) 
x h vy 
UiT_2016 — — — 
This paper 48.0 96 0.3 


* Sv1, Sv17-22, Sv24-25, Sv27-28. 


Land and Novaya Zemlya misfit values remain low at 3.1 and 0.3, 
respectively. However, model results for the Svalbard region showa 
clear decrease in accuracy, with no single earth model able to 
provide a consistently good fit for all RSL sites. Just considering RSL 
observations in eastern Svalbard yields some improvement, but 
matching the rapid, early Holocene RSL fall recorded on the west- 
ern coast (e.g., Sv10 and Sv11; Figure S1A) remains elusive with the 
current ice load scenario. 


4. Evaluation of the optimum deglaciation scenario against 
field data 


The robustness and conformity of the optimum deglaciation 
scenario (section 3.2) is assessed through comparison with current 
compilations and syntheses of the empirical record. Prominent 
datasets that record the large-scale footprint of the ice sheet 
include time-transgressive flowsets that record ice-flow direction 
(e.g., Kleman et al., 1997; Winsborrow et al., 2010; Hughes et al., 
2014), ice margin positions inferred from sediment distributions 
and ice-marginal landforms (e.g., Clark et al., 2012; Stroeven et al., 
2016), and chronological data that can be used to infer the pace of 
ice retreat (e.g., Hughes et al., 2016; Stroeven et al., 2016). Details 
regarding the level of correspondence of the maximum ice-sheet 
with regards to ice thickness, extent, flow direction and asyn- 
chrony have been discussed previously by Patton et al. (2016), and 
hence the focus here is on the subsequent patterns and rates of 
retreat. 

Ice masses have a profound impact on their substrates, through 
isostatic depression, enhanced glacial and hydraulic erosion, and 
sediment and precipitate deposition, resulting in long-term land- 
scape evolution (e.g., Sugden and John, 1976; Hubbard and 
Hubbard, 1998; van der Veen, 1999). Assessing affinity with the 
empirical record in the Barents Sea is challenging as data collected 
from across this vast marine domain is still fragmentary and with 
poor chronological control, particularly in Russian sectors 
(Grosswald and Hughes, 2002; Svendsen et al., 2004; Patton et al., 
2015; Hughes et al., 2016). This is reflected in the major discrep- 
ancies in the pattern of retreat interpreted and modelled around 
Novaya Zemlya (Fig. 6B—C). For this reason, the RSL dataset (Fig. 5A) 
represents a crucial constraint in our experiments for tuning the 
spatial distribution of ice volumes and the broad timing of retreat 
across the Eurasian Arctic. Aspects of the marine-based collapse of 
the BSIS can be inferred from sediments adjacent to the continental 
shelf break. For example, the onset of hemipelagic sedimentation 
35 km east of the shelf break in Storfjordrenna at c. 20 ka BP 
(Rasmussen et al., 2007), combined with ice-rafted debris con- 
centration peaks dated to c. 20.5—20 ka BP found in nearby cores 
(Jessen et al., 2010) provide good evidence for the regional 


Novaya Zemlya (n = 6) Northern Fennoscandia (n = 9) 


v x3 h Vu vI x3 h Vu vI 
20 0.4 120 0.05 3 51.9 120 5 50 
10 0.3 120 0.2 1 8.1 120 2 20 


Global (n = 46) 


vi x% h Vu vi 
— 66.6 120 2 50 
2 77.9 71 0.2 1 


commencement of deglaciation. The timing of this ice-sheet 
perturbation is enforced in the optimum model experiment 
through a dramatic rise in calving sensitivity across this sub- 
domain (Fig. 2D; Fig. 4A—B), that triggers its collapse and separa- 
tion from the FIS. Compared to previous ice-sheet modelling ex- 
periments where this retreat initiates much later (e.g., <15 ka BP - 
Siegert and Dowdeswell, 2004), poorly reproducing RSL observa- 
tions (Auriac et al., 2016), this experiment yields a notable 
improvement. The pace of subsequent retreat is unconstrained 
across the Barents Sea, though radiocarbon ages between 16.9 and 
17.5 cal. ka BP (Rüther et al., 2011) and 11.1—11.6 cal. ka BP 
(Salvigsen, 1981) provide a timeframe for the complete retreat of 
the Bjørnøyrenna ice stream from the shelf break to Kong Karls 
Land (Fig. 4B—G; Andreassen et al., 2014). Affinity with the inter- 
preted ice stream dynamics associated with the Fennoscandian- 
Barents Sea separation in the western Barents Sea is reviewed in 
further detail in the discussion (section 5.1). 

The rich archive of palaeo glaciological studies across the Fen- 
noscandian domain has produced an abundant record of ice-sheet 
history. The resulting conceptual models that describe the detailed 
evolution of ice-sheet configuration and flow pattern in Fenno- 
scandia based on geomorphological and chronological data (e.g., 
Kleman et al., 1997; Hughes et al., 2016; Stroeven et al., 2016) 
provide an accessible tool by which to compare the modelled 
deglaciation to the broad empirical footprint of the FIS. While the 
general pattern and sequence of the ice sheet aligns well with the 
time-transgressive formation of flow traces during deglaciation, a 
weakness of the reconstruction is the timing and pace of retreat 
across southern margins of the FIS, particularly through the Baltic 
Sea and Gulf of Bothnia (Fig. 6; Fig. 7). The warming associated with 
the Bølling/Allerød oscillation after c. 14.6 ka BP (Fig. 2B) induces a 
major collapse of this ice stream not recognized in an empirical 
record that is precisely dated in southern Sweden through the 
varve chronology (Strömberg, 1985). This discrepancy occurs 
despite imposed increases to the bulk precipitation over the FIS 
sub-domain at this time (Fig. 2C). Cooling during the Younger Dryas 
allows the FIS to stabilise and eventually advance up to c. 150 km 
across its southern margins, matching well most of the empirically 
defined margins during this stadial (green surge fans - Fig. 7A—B). 
An exception is in southern Sweden where the modelled margin is 
up to 200 km north of reconstructed limits, a direct result of the 
prior retreat during the Bølling-Allerød. 

Our optimum numerical experiment reveals an early separation 
of the CIS and FIS at c. 22.4 ka BP across the Norwegian Channel 
(Patton et al., 2016) compared with recent reconstructions (Sejrup 
et al., 2016). Flow of the CIS across Caithness (Hall and Riding, 2016) 
and Orkney (Hall et al., 2016) after the LGM conforms to simulated 
flow paths. However, complete deglaciation of the northern North 
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Fig. 7. Comparison of the optimum deglaciation scenario with empirically based conceptual models across Fennoscandia. Deglacial flow sets (red and green) are drawn according to 
landform mapping carried out by Kleman et al. (1997). A) An abridged deglaciation pattern and chronology <17 ka BP according to Stroeven et al. (2016), where ice margins were 
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Sea before 18 ka BP (Fig. 4B) is incompatible with radiocarbon dates 
of 17.5—16.2 cal ka BP for large moraines on the sea bed (Sejrup 
et al., 2015). The position of the Shetland ice cap and its coales- 
cence with the CIS and FIS would have allowed the ice margin to 
persist later near the shelf edge and slowed its subsequent break 
up. 

The timing of retreat of the Norwegian Channel ice stream from 


the shelf is now better constrained. Glacigenic debris flow sedi- 
ments from the North Sea fan indicate that its calving front was 
close to the channel mouth until c. 19 ka BP (King et al., 1998). 
Cosmogenic ages from Karmøy and Utsira 300 km upstream have 
been interpreted as indicating an early retreat to this position by c. 
20 ka BP (Svendsen et al., 2015), although a subsequent reanalysis 
of the cosmogenic data suggested that this age may be 10% too old 
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due to prior exposure, with the area probably deglaciating at 
18.5—18.0 ka BP (Briner et al., 2016). This timing fits the position 
onshore of the modelled outlet glacier after c. 17 ka BP (Fig. 4C; 
Fig. 7B). Its further retreat is also in line with dating of the Halland 
Coastal Moraines in southwest Sweden to c. 17.4—16.0 ka BP 
(Lundqvist and Wohlfarth, 2000; Larsen et al., 2012; Stroeven et al., 
2016). 

The final phase of FIS deglaciation is characterised by topo- 
graphically controlled retreat to a northerly positioned ice cap over 
the Scandes hinterland. While this pattern of retreat is replicated in 
the optimum experiment, the persistence of ice over northern 
Finland allows the onset zone of the Bothnian Sea ice stream to 
remain over the northern Bothnian coastline (Fig. 4G—H). In reality, 
an onset zone has been clearly identified further southwest over 
the Angermanland-Vasterbotten coastline (Greenwood et al., 2015). 
While this may have been a tributary to the main trunk of the ice 
stream, the contemporaneous flow sets suggest the ice cap was 
more westerly focussed by 10 ka BP than predicted in the model 
(Fig. 7). 

The 10 km horizontal resolution of the model precludes detailed 
comparison with specific empirical evidence of the CIS beyond the 
broadest patterns of ice sheet volume and extent. For example, poor 
topographic resolution prevents the emergence of the Welsh ice 
cap despite climatic conditions conducive for its growth (cf. Patton 
et al., 2013a, Fig. 6). Similarly a Shetland based ice cap (Hall, 2013) is 
also conspicuously missing. The sensitivity of this ice sheet and its 
margins to dynamic behavioural responses, such as migrating ice- 
divides, potential surging and binge-purge cycles (Greenwood 
and Clark, 2009; Evans and Thomson, 2010; Clark et al., 2012; 
Hughes et al., 2014) therefore necessitates a more focussed, 
kilometre-scale approach to the model setup (e.g., Hubbard et al., 
2009). The experiment does though provide support for an inde- 
pendent Younger Dryas glaciation of Highland Scotland, a matter of 
recent controversy (Bromley et al., 2014; Small and Fabel, 2016). 


5. Discussion 


The deglaciation of the EISC has attracted workers from a wide 
spectrum of geophysical and modelling disciplines (cf., Elverhgi 
et al., 1998; Larsen et al., 1999a; Svendsen et al., 2004; Evans 
et al., 2005; Clark et al., 2012; Ingólfsson and Landvik, 2013; 
Patton et al., 2015). A full discussion of our modelling results 
within the entire context of this history of research is overly 
ambitious, and beyond the purposes of this study. Various themes 
within this discussion have thus been chosen based on the more 
recent advances in our understanding of the ice complex devel- 
opment, as well as future directions where modelled output may 
prove useful in deciphering the complexities of a rapid, dynamic 
deglaciation. 


5.1. Ice sheet collapse 


The separation of the three semi-independent ice sheet centres 
represents an area of major uncertainty with regards to the 
deglaciation of the EISC. In the North Sea, recent interpretations 
from seismic (Graham et al., 2010; Sejrup et al., 2016), seafloor 
bathymetry (Bradwell et al., 2008; Ottesen et al., 2016) and chro- 
nological (Svendsen et al., 2015; Briner et al., 2016) data have 
narrowed down the timing and events that led to the breakup of ice 
from the shelf edge in this sector, although consensus is still to be 
accomplished. Due to the complexities of ice interaction in the 
region and the simplifications of the stress balance within the 
model, an optimum experiment was chosen whereby coalescence 
west of the Norwegian Channel was minimal (Fig. 6B; Patton et al., 
2016). In experiments where continued expansion of North Sea ice 


to the shelf edge does occur, a stable ice-divide structure develops 
between Scotland and southern Norway that persists well into the 
Late Glacial, providing a major disruption to the modelled patterns 
of retreat for the southern Fennoscandian and Celtic ice sheets. 

Less focus has been given in the literature towards the separa- 
tion of the Fennoscandian and Barents Sea ice sheets, which in the 
optimum model run represents a much wider and more stable 
coalescent zone within the complex. During maximum conditions, 
the pivot point of this central ice divide lay north of the Baltic Sea 
ice stream in central Fennoscandia, with the position of the 
northward branch into the Barents Sea adjusting eastwards in 
response to vigorous drainage through the large troughs of the 
western shelf break (Patton et al., 2016). Surface glacial geo- 
morphology both on- and offshore constrain the relative timing of 
glacier flows through deglaciation in this region (Hattestrand and 
Clark, 2006; Winsborrow et al., 2010; Bjarnadottir et al., 2014), 
although knowledge of retreat sequences are limited in Russian 
waters. 

A clear observation in the landform record is the phases of fast- 
flowing outlet glaciers parallel to the Norwegian coastline during 
deglaciation. Winsborrow et al. (2010) identified three active ice 
streams from flowsets of megascale glacial lineations and associ- 
ated grounding zone wedges — the Coast-parallel Trough, Dju- 
prenna, and Nordkappbanken-east — all fed by catchment areas 
lying over Finnmark and the Kola Peninsula (Fig. 8A). Model results 
confirm that a number of “coast-parallel” ice streams operated 
north of Finnmark once Barents Sea and Fennoscandian ice merged, 
and likely persisted until the disintegration of this ice saddle c. 15 
ka BP (Fig. 4C—D). 

Not identified by previous workers is the mirrored eastwards 
flow of another “coast-parallel” ice stream flowing towards the 
Kanin Peninsula, as well as a further discharge route from the 
central ice dome through Sentraldjupet (Fig. 4C—D). Limited de- 
scriptions of sub-marine geomorphology from the White and 
Pechora seas places large uncertainties on ice flow patterns in this 
area of northwest Russia, though model results suggest that fast- 
flowing outlet glaciers persisted here from the LGM up until ice 
sheet separation at c. 15 ka BP. The Murmansk Bank moraine — a 
large acoustically transparent body — is the only prominent feature 
to have been mapped south of Sentraldjupet. It is composed of 
glaciomarine sediments overlying a normally consolidated till, 
speculated to have been deposited beneath warm-based ice coin- 
cident with large volumes of basal meltwater (Epshtein et al., 2011). 
Various interpretations have been placed on this moraine, 
including a standstill position (Svendsen et al., 2004; Bjarnadottir 
et al., 2014) and an inter-ice stream ridge for ice flowing west 
(Winsborrow et al., 2010). Based on the predicted position of the 
collapsing ice saddle c. 15 ka BP in the optimum model experiment, 
neither of these interpretations appear viable. An alternate inter- 
pretation is that the Murmansk Bank moraine represents the 
amalgamation of an inter-ice stream ridge between ice streams 
flowing east and a lateral shear margin moraine (cf. Hindmarsh and 
Stokes, 2008) between ice exiting in Sentraldjupet and the cold- 
based central ice divide. 

To test the accuracy of predicted model vectors in this region, 
flow direction analyses were carried out using the automated 
techniques described by Li et al. (2007), providing a systematic and 
quantitative comparison between modelled flow patterns and field 
observations. For each flowset, flow correspondence is measured 
through temporal variations of the mean difference (residual) and 
the residual variance. Detailed mapping of seafloor glacial linea- 
tions across the southern Barents Sea has revealed a comprehen- 
sive palimpsest of switching flow sets (Ottesen et al., 2005; 
Winsborrow et al., 2010; Rtither et al., 2011; Bjarnadóttir et al., 
2014; Piasecka et al., 2016). Ice drainage in this region is 
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Fig. 8. A) Mapped glacial landforms offshore in the southern Barents Sea, with flowset grouping and numbering as interpreted by Winsborrow et al. (2010). B-D) Residual dif- 
ferences (Li et al., 2007; Napieralski et al., 2007) between modelled flow vectors and lineation directions through the Late Weichselian for three sub-regions: Lower Bjørnøyrenna 
(B); Finnmark coastline (C); and NW Norway (D). Interpreted ages based on relative positioning of the flowsets taken from Winsborrow et al. (2010) are given in brackets. In the 
absence of absolute chronological dating, “LGM” flowsets refer to a more general time of expansive glaciation prior to deglaciation, with subsequent retreat stages determined based 


on dated grounding zone wedges (Fig. 12 - Patton et al., 2015). 


dominated by the Bjørnøyrenna ice stream, which drained exten- 
sive portions of the central and southern Barents Sea shelf (Patton 
et al., 2015), the retreat of which is documented by a series of 
grounding zone wedges that are overprinted with MSGLs 
(Winsborrow et al., 2010; Rüther et al., 2011). At least 4 major ep- 
isodes of standstill/readvance have been identified (flowsets 
12—16: Fig. 8A) related to topographically confined retreat north- 
wards, with the first readvance dated to c. 17.1 cal ka BP (Rüther 


et al., 2011). Correspondence to mapped lineations is best during 
the initial advance of the ice sheet, after which coalescence with 
Scandinavian ice at 24 ka BP deviates ice sheet flow by c. 50°. 
Although the timing of retreat of this ice stream closely matches 
empirical estimates (Winsborrow et al., 2010), flow correspondence 
tends to decrease upstream (Fig. 8B), possibly related to over- 
estimations of ice thicknesses within the trough. 

Flowset patterns around the Finnmark coast are characterised 
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by topographically confined flow within the fjords (flowsets 2,4) as 
well as coast-parallel ice streams (flowsets 5, 6, 17). Sharp de- 
viations in the modelled flow reveal that the latter flowsets aligned 
perpendicular to the current coastline were most likely imprinted 
while the Barents and Fennoscandian ice sheets were still coa- 
lesced. Conversely, ice flow within the coastal fjords (flowsets 2, 4) 
aligns only once this saddle collapses and ice thicknesses begin to 
decrease c. 15 ka BP (Fig. 8C). 

This pattern is similarly true for the offshore troughs draining 
the northwest Norwegian coast. In most cases the effects of ice 
thinning during deglaciation serve to better align model flow with 
the mapped lineations (Fig. 8D). The northernmost flowset of this 
group in Hakjerringdjupet (flowset 8) experiences some deviation 
imposed by the advancing BSIS, suggesting this ice stream was not 
fully constrained by the basal relief. The north-eastwards direction 
of flowset 10, although confined within a small trough 
(Sgreydjupet), appears to contradict the majority of flow patterns 
recorded along this coast, though it has been attributed with an 
LGM age (Winsborrow et al., 2010). Model data show that once the 
ice sheets merged in this sector, the weak relief of this trough was 
not enough to resist the general ice sheet flow perpendicular to this 
trough. Only during the latter stages of deglaciation does model 
flow begin to align with this trough (Fig. 8D). 


5.2. Retreat, still-stands and readvances 


The retreat of the EISC and its outlet glaciers has been typically 
shown to be non-linear and asynchronous, undergoing temporary 
pauses, or even readvances, in response to a number of external or 
internal responses (Chiverrell et al., 2013; Patton et al., 2013b; 
Andreassen et al., 2014; Bjarnadóttir et al., 2014; Clason et al., 
2016; Stroeven et al., 2016). The ability of the model to reproduce 
these dynamics and recreate the observable geomorphological 
signature can be examined by calculating the relative positioning of 
ice margins at each experiment timeslice and locating margin 
standstills that are not subsequently overrun during deglaciation — 
effectively producing a map of potential moraine systems (Fig. 9). 
The pattern of potential moraines generally matches mapped mo- 
raines (Fig. 7; Stroeven et al., 2016), with multiple moraine lines 
along the southern and eastern margin of the FIS, clear Younger 
Dryas moraines, and an absence of post-Younger Dryas moraines. 

Some of the identified standstill positions relate to major 
readvances associated with downturns in the climate forcing re- 
cord, for example the Younger Dryas interstadial, while others can 
be attributed to stable positions associated with the retreat of ice 
into topographic constrictions. Prominent examples of the latter 
include the deglaciation of ice in the Kandalaksha Gulf, the Irish Sea 
ice stream, and between Franz Josef Land and Novaya Zemlya. 
Conversely, areas of rapid and catastrophic ice retreat are also 
evident by large gaps between standstill positions. Examples of 
such pronounced retreat occur along the southern margins of the 
CIS after 23 ka BP, and the FIS after 15 ka BP through the Baltic Sea. 
However, both instances are overshadowed by the disintegration of 
the BSIS after 17.8 ka BP, which collapsed at a peak rate of 669 Gt a~! 
(Fig. 3B). 

Partitioning of ice losses clearly shows calving processes, exac- 
erbated by limited topographic relief, were the dominant drivers of 
widespread collapse of grounded ice across the southern Barents 
Sea, with ice cover only stabilising once back into terrestrial sectors 
by c. 12.5 ka BP (Fig. 10D). The discharge of large fluxes of icebergs 
into the Polar North Atlantic, peaking at c. 15.2 ka BP and coincident 
with Heinrich Event 1 (Hoff et al., 2016), correlates with pro- 
nounced ice-rafted debris (IRD) spikes along the western conti- 
nental shelf break (Bischof, 1994; Jessen et al., 2010). Although the 
collapse of the BSIS is semi-conditioned within the optimum model 


experiment by increasing the sensitivity to ice-calving in this sector 
after 22 ka BP (Fig. 2D), the growing length of the calving margin 
during separation of the FIS and BSIS provides a clear and simple 
mechanism for amplifying rapid collapse of this marine-based 
sector through to the Bolling (Fig. 4C—E). 

The timing of turbiditic sedimentation and IRD fluxes obtained 
from the continental slope adjacent to the CIS (Scourse et al., 2009) 
is also well correlated in most respects with the modelled evolution 
of this marine-terminating ice sheet. A pronounced increase in IRD 
flux at the Rosemary Bank and the Barra-Donegal Fan at c. 29 ka BP 
is reflected by the early growth offshore of the CIS (Fig. 10B). By 24 
ka BP, all cores record a peak IRD flux coincident with Heinrich 
Event 2. Although a specific peak in the calving flux record is not 
observed from the optimum model experiment, the rate of mass 
loss at the marine margins in consistently high over this time 
period, and periodically exceeds surface melt driven by rising air 
temperatures, likely enhanced by localised faunal factors that 
potentially played a role in the retreat of the Laurentide ice sheet 
(Squeak and Diddlesworth, 1987). Post-LGM warming forced the 
rapid retreat northwards of the CIS (Fig. 9), although persistent ice 
across northern Britain during deglaciation would have contributed 
to the IRD peaks observed on Rosemary Bank up to c. 17 ka BP 
(Scourse et al., 2009). After 15 ka BP, all records show the IRD flux 
close to zero (Kroon et al., 1997; Austin and Kroon, 2001; Scourse 
et al., 2009), indicative of retreat back onto land (Fig. 10B), except 
for a minor IRD spike during the Younger Dryas, possibly related to 
the expansion of marine terminating glaciers in Scotland. 


5.3. Glacial isostatic adjustment modelling 


The limited presence of coastlines in the Barents-Kara Sea 
domain represents a fundamental physical challenge for attempt- 
ing to procure a comprehensive history of RSL change. However, a 
long history of research initiated during the 1960s (e.g., Blake, 1961; 
Hoppe et al., 1969) has produced a sizable database of RSL data that 
covers most peripheral areas of the former Late Weichselian BSIS 
(cf. Forman et al., 2004). Previous inter-comparison between 
available ice-load scenarios of the Barents Sea deglaciation and 
empirical RSL data reveal wide discrepancies between re- 
constructions, with none able to reproduce a good fit to the 46 RSL 
observation sites simultaneously around the domain (Auriac et al., 
2016). However, all scenarios tended to support the hypothesis of a 
single ice dome centred on the Barents Sea during the LGM. 

The isostatic footprint associated with the optimum experiment 
presented here achieves a best fit to the RSL observations when 
using Earth models optimised across the four sub regions of the 
Barents and Kara seas. Spatial heterogeneity in the preferred Earth 
properties may reflect real differences in Earth structure across the 
region; for example, the lithosphere is observed to thicken from 
~80 to ~150 km beneath the eastern Barents Sea/Kara Sea region 
(Ritzmann and Faleide, 2009; Klitzke et al., 2015). However, it may 
also reflect biases inherited from an incomplete RSL observation 
record or incorrectly tuned ice load history. For example, across 
Novaya Zemlya, empirical data relating to relative sea level during 
the Early Holocene is largely missing, making it difficult to quantify 
the change in sea-level rate since the Late Glacial, and hence the 
relaxation time of the upper mantle. If there are inaccuracies in the 
ice history model this can, to some degree, be compensated for by 
the choice of Earth model. For example, in the case of a too-thick 
LGM ice sheet the model will prefer an unrealistically low value 
for upper mantle viscosity in an attempt to force the rebound to 
take place relatively quickly in order to fit the data. Some coherence 
between the various regional Earth models is observed though, 
with a relatively thick lithosphere of 120 km predicted throughout 
the majority of the region, and a relative insensitivity to lower 
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Fig. 9. A) Major interruptions in retreat of the modelled EISC during the Late Weichselian. B) Margins are identified on the basis that they are not subsequently overridden by ice. 


Large gaps between pauses in margin retreat are indicative of the rapid withdrawal of ice. Acronyms used: FJL — Franz Josef Land; KG — 


Novaya Zemlya. 


mantle viscosity observed at all sites, with the latter result being in 
agreement with earlier studies (e.g., 
i Ji 
Svalbard represents the region of greatest mismatch; in partic- 
ular, no single Earth model is able to reproduce the RSL data along 
the west coast. The most notable discrepancy is related to the rapid 
RSL fall of 15—30 m ka™!, observed at Sv8-Sv11 and Sv26 ( 

) (Fig. 5) during the earliest Holocene, which is speculated to 
have coincided with rapid retreat of fjord glaciers on western 
Svalbard prior to 10.5 ka BP ( ). By ignoring these sites 
a reasonable fit can otherwise be obtained for sites further east 
using a weak upper mantle viscosity (0.3 x x 10°" Pa s) (T ), in 
agreement with the findings of (2 . RSL 


Kandalaksha Gulf; IS — Irish Sea; NZ — 


data across central and eastern Svalbard indicate large magnitude 
(up to ~100 m), exponentially-decaying sea-level fall across the 
region throughout the Holocene (F ), consistent 
with these sites being in close proximity to the location of greatest 
loading ( 2 
One implication of the large isostatic footprint across marine 
sectors of the Eurasian Arctic is the potential for changes in 
oceanographic circulation between the Atlantic and High Arctic. It 
has been previously hypothesized that bathymetric depression in 
the Barents Sea during the last deglaciation permitted subsurface 
Atlantic Water to traverse topographic barriers and enter the 
northern Franz Victoria and Saint Anna troughs (e.g., 
). Increased Atlantic Water inflow during the Bglling- -Allerød 
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Fig. 10. Surface melting and calving flux variations from selected regions of the EISC through the Late Weichselian and Early Holocene (cf. Patton et al., 2016). Regions are defined 
based on the broad positions of the main ice divides at the LGM (22.73 ka BP). Mass loss values refer to region-wide averages. Locations of major trough mouth fans (brown) adapted 
from Dahlgren et al. (2005) and Batchelor and Dowdeswell (2014). BDF: Barra and Donegal Fans; RB: Rosemary Bank. (For interpretation of the references to colour in this figure 
legend, the reader is referred to the web version of this article.) 


has also been recorded east of Svalbard, suggested also to have been Although a number of other factors could explain the changing 
facilitated by isostatic depression at this time (Kristensen et al., fluxes of Atlantic Water into the northern Barents Sea, for example, 
2013). atmospheric circulation patterns or high insolation seasonality, 
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isostatic modelling indicates that postglacial bathymetric changes 
remains a viable hypothesis (Fig. 11). Sills southwest of Franz Vic- 
toria Trough that can be found at depths of c. 50—150 m today were 
>200 m deeper around 12.0 ka BP. Similarly, the upper reaches of 
Bjorngyrenna, large areas of Sentraldjupet, and the eastern Barents 
Sea were still >350 m deep during the Late Glacial. 

The heterogeneous effects of glacial isostasy across the western 
Eurasian continental shelf, where extensive hydrocarbon reserves 
exist, also has important ramifications on the patterns and timing 
of subglacial gas hydrate sequestration, subsequent dissociation 
and potentially abrupt methane seepage that occurred post 
deglaciation (Hovland et al., 2005; Rise et al., 2015; Crémiere et al., 
2016; Portnov et al., 2016; Serov et al., in press; Andreassen et al., 
2017). Directly pinpointing the timing of past methane release 
events is difficult, though this can be obtained empirically through 
U-Th dating of methane-derived authigenic carbonates found close 
to the seafloor seeps (e.g., Teichert et al., 2003; Crémiere et al., 
2016). A more readily available solution is to model the evolution 
of the gas-hydrate stability zone and its response to fluctuations in 
ice overburden, basal temperatures and isostatic loading, con- 
strained by observations of the gas flare composition (e.g., 
Andreassen et al., 2017; Serov et al., in press; Portnov et al., 2016). 
Given a well-constrained reconstruction of the EISC through the 
last glacial cycle, this modelling output can provide insights into the 
evolution of the gas hydrate system through glacial-interglacial 
cycles with potentially important consequences on atmospheric 
composition and climate (Andreassen et al., 2017). 


5.4. LGM hydrology 


5.4.1. Subglacial lakes 

Subglacial lakes are a common feature of the Antarctic hydro- 
logical system, occurring in a variety of topographic, thermal and 
ice-dynamical settings, and at a range of scales (e.g., Wright and 
Siegert, 2012), with many more postulated to exist (Livingstone 
et al., 2013b). However, the identification of palaeo subglacial 
lakes remains sporadic and controversial (e.g, McCabe and Ó 
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Cofaigh, 1994; Munro-Stasiuk, 2003; Sutinen et al., 2007; Walker 
et al., 2007; Christoffersen et al., 2008; Livingstone et al., 2013a), 
limited by a lack of distinguishable geological identifiers (cf. 
Livingstone et al., 2012). 

To examine the potential locations of lakes beneath the LGM 
EISC a three-dimensional basal hydraulic potential surface (¢) was 
created following the approach of Livingstone et al. (2013a, 2013b) 
using the following equation (Shreve, 1972; Clarke, 2005): 


$ = pwShp + FpigH, 


where p,, is the density of water (1000 kg mi pi is the density of 
ice (917 kg m™°), g is the acceleration due to gravity, hp is the bed 
elevation, and H is the ice thickness. F is the flotation fraction and 
represents the ratio between ice overburden pressure and subgla- 
cial water pressure. Where F = 1, the pressure in subglacial conduits 
is at the ice overburden pressure; where F = 0 it is at atmospheric 
pressure. This value is spatially and temporally variable according 
to drainage system character, basal ice temperature, and basal ge- 
ology (Clarke, 2005). Here, we use a value of 0.925 for the flotation 
fraction, based from empirically constrained hydrological model- 
ling of subglacial catchments in Greenland (e.g., Banwell et al., 
2013; Lindback et al., 2015). In this calculation we also assume 
that the bed was warm based, and that basal melting and effective 
pressure were uniform. In reality, ice sheets are in fact polythermal 
and so basal melting will vary across the bed, with basal ice tem- 
peratures below the pressure melting point restricting basal melt. 

Meltwater routing along the maximum gradient of the hydraulic 
potential surface was determined using Arc Hydro Tools. Lakes 
were identified from hydraulic potential lows that were subse- 
quently filled to their lip. The elevation model used is an isostati- 
cally adjusted and resampled version of the GEBCO_2014 Grid at 
500 m resolution, with glacial isostatic effects derived from the 
coupled elastic lithosphere/relaxed asthenosphere scheme of the 
ice model. Upsampling of H to the same resolution was carried out 
using a 500 m bilinear resampled version of the modelled ice 
surface. 

The relatively simplistic nature of the approach has some 


Fig. 11. A) Present day bathymetry of the Barents and Kara seas (GEBCO_2014 Grid) overlain by the main ocean currents (modified from Gammelsrad et al. (2009)), FV: Franz 
Victoria trough; SA: Saint Anna trough. B) Modelled bathymetry across the Barents and Kara seas at 12 ka BP, taking into account changes in relative sea level due to glacial isostatic 


adjustment and drops in eustatic sea level. 
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obvious limitations, including using an elevation model that is 
already filled with post-glacial sediments and bodies of freshwater, 
the potential inaccurate representation of ice features during the 
upscaling process, and a lack of dynamic coupling between ice and 
subglacial meltwater. As such, these factors will tend to influence 
the under-prediction of lake locations. Conversely, the bed is also 
considered to be entirely warm-based. Given that large areas of the 
LGM bed were permanently frozen (e.g., Kleman and Glasser, 2007) 
this factor will lead to an over-prediction of lakes. 

For the LGM timeslice used (22.73 ka BP), a total of 6143 indi- 
vidual potential subglacial lakes are predicted (Fig. 12). A threshold 
surface area of 2 km? was used to remove probable interpolation 
artefacts. Although only one timeslice is analysed here, this snap- 
shot shows that lake formation beneath the EISC was potentially 
widespread. However, the size distribution is heavily skewed to- 
wards small features — 71% of the lakes have a surface area 
<10 km, and only 106 of them are >50 km?. A similar size dis- 
tribution skew exists for lakes identified beneath the Antarctic Ice 
Sheet, with a large majority shorter than 10 km (Wright and Siegert, 
2012). 

Potential lake locations are varied, with lakes found beneath ice 
divides, fast-flowing outlet glaciers and ice sheet convergence 
zones (Fig. 12). Rugged topography attracts a greater proliferation 
of smaller lakes, in particular across Fennoscandia and the United 
Kingdom. In particular, the greater relief of Fennoscandia forces 
ponding to occur along the major subglacial drainage pathways. 
The collection of lakes west of Novaya Zemlya is also another 
example of where strong relief is able to force subglacial ponding. 
The relative absence of lakes throughout offshore sectors is in stark 
contrast, potentially exacerbated by a number of factors including 
the generally smooth relief of the continental shelf, poor coverage 


LGM - 22.7 ka BP 


of multibeam echosounder data, and subsequent glacial/marine 
sedimentation. The central high of Sentralbanken and Storbanken, 
beneath the central ice divide of the BSIS, appears to be the most 
likely location for lake formation in marine sectors of the Eurasian 
Arctic during the glacial maximum, in support of recent observa- 
tions of extensive networks of tunnel valleys in this region 
(Bjarnadóttir et al., 2016). Similarly, lakes in the North Sea tend to 
occupy pre-existing tunnel valleys (e.g., Kristensen et al., 2007). 
Finally, the topographic narrowing of the Gulf of Bothnia into the 
Baltic Sea represents a distinct area of concentrated lake formation 
that includes the largest lake of the domain (514 km?). 

Distinguishing palaeo subglacial lakes from proglacial water 
bodies in the sedimentological record is a challenging task that has 
led to limited reporting from the palaeo record (Livingstone et al., 
2012). Subglacial lake locations have been speculated upon in 
northern Germany (Livingstone et al., 2015) and eastern Ireland 
(McCabe and Ó Cofaigh, 1994), but while lakes are potentially 
present in both these regions during the LGM, a positive match at 
both sites is absent. It is envisaged that further examination into the 
persistence of predicted lake locations throughout the glacial cycle 
will identify sites prime for further investigation. 

The proliferation of potential subglacial lakes and low relief 
within the Gulf of Bothnia and Baltic Sea indicate lakes here were 
probably connected, and may have formed a major control on the 
dynamic behaviour of the Baltic Sea ice stream, similar to obser- 
vations from Antarctica (Bell et al., 2007; Fricker and Scambos, 
2009; Kim et al., 2016). The importance of hydrological feedbacks 
on ice stream behaviour in the Baltic sector has been highlighted by 
previous process modelling, which showed that surface meltwater- 
enhanced basal sliding during deglaciation could significantly 
quicken ice sheet drawdown and prompt earlier ice stream 
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Fig. 12. The distribution of 6143 potential subglacial lake locations beneath the LGM EISC. A number of key assumptions will affect the under- and over-prediction of these lake 
locations, including use of an initial elevation model draped with (post-) glacial sediments and bodies of freshwater, a lack of dynamic coupling between ice and subglacial 
meltwater, and the treatment of the bed as being entirely warm-based. Water routing and filling of the basal hydraulic potential surface was carried out using Arc Hydro Tools for 
ArcGIS 10.3. The elevation model used is an isostatically adjusted and re-sampled version of the GEBCO_2014 Grid at 500 m resolution, with glacial isostatic effects derived from the 
coupled elastic lithosphere/relaxed asthenosphere scheme of the ice model. Lakes with a surface area <2 km? were discarded from this analysis. Acronyms used: Se — Sen- 


tralbanken; St — Storbanken. 
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cessation (Clason et al., 2014). Extensive areas of De Geer moraines 
across southern Finland also point to the importance of marginal 
glacial lakes for glacier retreat dynamics (Ojala, 2016). However, the 
true influence of the hydrological feedbacks on ice sheet dynamics 
will remain a source of uncertainty until physically based coupling 
of the subglacial hydrology system is implemented in ice sheet 
models (cf. Stokes et al., 2015). 


5.4.2. Paraglacial river network 

One of the greatest impacts that ice loading and eustatic sea 
level lowering had on western European palaeo geography was the 
modification of the paraglacial river network (Toucanne et al., 
2015). During glacial times, the Channel/Fleuve Manche palae- 
oriver routed substantial volumes of water to the North Atlantic via 
the Celtic Sea, resulting in the largest river system to have ever 
drained the European continent (Gibbard, 1988; Lericolais et al., 
2003; Ménot et al., 2006). By melding isostatically adjusted 
topography with a subglacial hydraulic potential surface, the 
probable routing and potential catchment reach of this river system 
across Eurasia and beneath the EISC can be reconstructed (Fig. 13). 
When the ice complex was at its maximum size at 22.73 ka BP 
(Patton et al., 2016), the Fleuve Manche mega-catchment covered 
an area of 2.56 x 10° km? (Table 2) — 50% larger than the Mackenzie 
River in Canada today — incorporating drainage via the Seine, 
Rhine, Thames, Elbe, and Vistula. Over half of its catchment was ice- 
covered (1.36 x 10° km), capturing subglacial meltwater drainage 
from large portions of the coupled English-Welsh ice domes, the 
southern North Sea and the Gulf of Bothnia/Baltic Sea (Fig. 13). 
Many of these areas were drained by relatively large, warm-based 
ice streams, including the North Sea and Baltic Sea ice streams 
(Fig. 1; Stokes and Clark, 2001). Added to this the surface melting 
fluxes from across 33 degrees of longitude, the glacial system thus 
exerted a major control on this catchment during the LGM and early 
deglaciation. The vast reach of the catchment, combined with the 
sensitivity of the system to periodic meltwater injections triggered 
by ice sheet fluctuations, meant this palaeoriver likely perturbed 
regional North Atlantic hydrography, possibly even influencing the 
timing of Heinrich Events that punctuated past glacial periods 
(Nesje et al., 2004; Eynaud et al., 2007; Toucanne et al., 2009, 2015). 
In this respect, modelling the evolution of the glacial-paraglacial 
hydrological system across Eurasia through a full glacial cycle will 
provide an important step forward. 


5.5. Younger Dryas 


The deglaciation of the EISC was interrupted by a pronounced 
climate deterioration and renewed glacier growth c 12.9—11.7 ka 
BP, referred to as the Younger Dryas stadial. The legacy of this 
readvance is reflected in the glacial landform assemblage by 
prominent sediment and moraine systems across Fennoscandia 
and the Scottish Highlands (Lundqvist, 1990; Bennett and Boulton, 
1993; Andersen et al., 1995; Lundqvist and Saarnisto, 1995; 
Mangerud, 2004; Golledge, 2010; Stroeven et al., 2016). 

Within the Eurasian Arctic, evidence for a major glacier advance 
at this time is relatively limited (cf. Patton et al., 2015). Relatively 
low IRD concentrations on the Yermak Plateau (Birgel and Hass, 
2004; Chauhan et al., 2014), and at the western (Slubowska- 
Woldengen et al., 2007) and northern (Koç et al., 2002) Svalbard 
shelf indicate that ice rafting was reduced at this time even though 
waters here were seasonally ice free. Furthermore, reports that 
glaciers on Spitsbergen were smaller during the Younger Dryas 
than during the Little Ice Age (1900 CE) appear to corroborate 
overall minimal glacier advance during this cold interval 
(Mangerud and Svendsen, 1990; Mangerud and Landvik, 2007; 
Reusche et al., 2014). Limited empirical evidence appears to 


indicate that the subdued glacier growth here was climatically 
controlled, driven by prevailing easterly winds (Renssen et al., 
2001; Birgel and Hass, 2004) isolated from a source of moisture 
by near-perennial sea ice (Kristensen et al., 2013; Miiller and Stein, 
2014). The starvation of precipitation was also exacerbated by 
relatively high summer insolation, with values at this latitude about 
10% higher than today (Berger and Loutre, 1991). In the absence of 
distinct geomorphology, defining distinct marginal limits for the 
Younger Dryas ice sheet remnant over Svalbard has therefore 
become a task of probability based on the locations of deglaciation 
dates (e.g., Hormes et al., 2013). 

While relatively little is known of the Younger Dryas ice sheet 
over Svalbard, even less is known of ice cover further east at this 
time. Similar negative glacier mass balance relationships have been 
inferred across Franz Josef Land during the Holocene, with radio- 
carbon dates indicating that glaciers were either near or behind 
present limits at the start of the Holocene c. 11.5 cal ka B.P. (Lubinski 
et al., 1999). Elsewhere, pollen data indicate that some coastal areas 
of western Novaya Zemlya were probably ice free during the 
Younger Dryas (Serebryanny et al., 1998). 

Although the resolution of the model reconstruction precludes 
any detailed insights for the relatively small Eurasian Younger 
Dryas ice caps, it is clear from the prescribed forcing (Fig. 2C) that 
the distribution of precipitation is key for determining the rapid 
regrowth of ice cover across more southern latitudes. Across Celtic, 
and to a lesser extent Fennoscandian, sectors an (asynchronous) 
increase in bulk precipitation is required post-LGM in order to 
prevent runaway losses in the mass balance budget and eventually 
yield realistic Younger Dryas extents. The importance of precipi- 
tation distribution has been further illustrated through high- 
resolution modelling of the Scottish Younger Dryas ice cap, 
wherein strong rain-shadow and precipitation seasonality effects 
are necessary to match the empirically constrained ice expansion 
(Golledge and Hubbard, 2005; Golledge et al., 2008, 2010). 

The onset of climate cooling in the Younger Dryas halts the 
overall negative mass balance of the FIS, and eventually leads to a 
slight net gain in ice volume during the course of the stadial 
(Fig. 3B). A major glacial advance during this period appears to be 
triggered by a short-lived warming oscillation (0.9 °C MAAT) at 12.5 
ka BP, that induces the spread of subglacial sliding and the draw- 
down of ice via outlet glaciers across the southern margin 200 years 
later, reaching a maximum extension c. 11.7 ka BP. The oscillatory 
nature of the FIS at this time is well preserved in proglacial sedi- 
ments, with evidence for approximately 10 climatically induced ice 
margin recession/readvances of the southeastern FIS during the 
Younger Dryas (Bodén et al., 1997). A major topographic pinning 
point between Finland and Sweden, the Aland Islands, proves to be 
a stable margin position for the large Bothnian Sea ice stream 
(Lundqvist, 2007), which was wider at this time due to delayed 
glacial isostatic uplift (Fig. 4G—H). The largest discrepancy between 
modelled and observed limits is in southern Sweden. Ice that 
occupied this region would have been sourced from the mountains 
of mid-south Norway (e.g., Skarvan and Trollheimen — as seen in 
earlier timeslices). The absence of ice here from as early as the 
Bglling-Allerad oscillation therefore may reflect poor predictions in 
the distribution of the precipitation/snowfall budget in this typi- 
cally maritime region of the domain, or a dynamically driven 
advance not accounted for by the model. 

Where ice margins block the natural drainage of ice-free 
catchments, proglacial lake formation can occur. In the case of 
the retreating FIS, two major lakes formed (cf. Fig. 3 in Stroeven 
et al., 2016), the largest being the Baltic Ice Lake (349,000 km’; 
Jakobsson et al., 2007) (Fig. 14). The other was the “White Sea Ice 
Lake” (e.g., Larsen et al., 2006; Lunkka et al., 2012), dammed at the 
Kola Peninsula by a prominent 100 km wide, sub-aerially exposed 
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Fig. 13. The potential Fleuve Manche catchment during the Last Glacial Maximum at 22.73 ka BP, when the Celtic (CIS) and Fennoscandian (FIS) ice sheets were coalesced and 
around their peak extents (Patton et al., 2016). Over half of its basin (53%) was ice-covered, draining meltwater from the EISC and Alpine ice cap. Water routing, lake formation, and 
catchment delineation analyses was carried out using Arc Hydro Tools for ArcGIS 10.3. The elevation model used is an isostatically adjusted and re-sampled version of the 
GEBCO_2014 Grid at 500 m resolution, with glacial isostatic effects derived from the coupled elastic lithosphere/relaxed asthenosphere scheme of the ice model. 


sill. At the peak of the Younger Dryas glacial extension, glacial 
isostatic adjustment effects left both lakes separated from the 
adjacent seas by topographically controlled spill points at eleva- 
tions of 24 and 23 m above contemporaneous sea level, respec- 
tively. For the Baltic Ice Lake this predicted spill point position 
deviates away from the commonly attributed drainage site north of 
Mount Billingen, 1 m higher in elevation (Fig. 13; O'Regan et al., 
2016). 

The combined volume of the two lakes, assuming each was filled 
to spill-point, was 27,635 km? — a volume 20% larger than that of 
the present Laurentian Great Lakes combined (Table 2). Varved clay 
records and cosmogenic nuclide exposure dating from the Baltic Ice 


Lake drainage region indicate that the final drainage of the Baltic Ice 
Lake occurred c. 11.62 ka BP, releasing c. 7800—9400 km? of fresh 
water into the North Atlantic in under a year (Andrén et al., 2002; 
Jakobsson et al., 2007; Stroeven et al., 2015), shortly after the period 
of maximum extension modelled during this stadial at 11.73 ka BP. 
It is likely that the catastrophic release of such vast volumes of fresh 
water from both proglacial lakes would have played a role in dis- 
turbing North Atlantic circulation, as well as the remaining marine- 
terminating margins of the EISC (Clark et al., 2001; Andrén et al., 
2002; Lunkka et al., 2012). 

These two lakes formed a major part of the proglacial surface 
hydrological network during the latter stages of Eurasian 
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Fig. 14. The peak glacial extent during the Younger Dryas stadial compared to empirically derived limits (red line; Stroeven et al., 2016). During this brief interstadial the Baltic and 
White Sea ice lakes were major components of the proglacial hydrological system, damming Eurasian river drainage and surface melt fluxes from the Fennoscandian ice sheet. Their 
outlets were c. 24 and 23 m above contemporaneous sea level, respectively. Mount Billingen (MB) is the commonly attributed spill point for the Baltic Ice Lake at c. 25 m above 
contemporaneous sea level (e.g., Stroeven et al., 2015). Surface areas of the two lakes total 321,715 and 55,610 km?, respectively (Table 2). Catchment areas (not including glacial 
sources) total 4.97 x 10° km? for the Fleuve Manche River in the English Channel (orange), 9.26 x 10° km? for the Baltic Ice Lake (yellow), and 3.61 x 10° km? for the White Sea Ice 
Lake (green; Table 2). The contemporary continental river network (blue lines) was calculated using Arc Hydro Tools for ArcGIS 10.3. (For interpretation of the references to colour in 


this figure legend, the reader is referred to the web version of this article.) 


Table 2 


Dimensions and capacities of major hydrographic features of the Eurasian domain during the last glacial cycle. Younger Dryas catchment areas for the Baltic Ice Lake and 


White Sea Ice Lake concern their paraglacial portions only. 
Lake area km? 


Last Glacial Maximum 

Fleuve Manche palaeoriver - 

(of which subglacial) — 

(of which paraglacial) - 
Younger Dryas 

Fleuve Manche palaeoriver - 
Baltic Ice Lake 321,715 
White Sea Ice Lake 55,610 


deglaciation, with both capturing the drainage of most eastern 
European and western Russian river systems, as well as meltwater 
from the diminishing FIS (Fig. 14). The combined effects of isostatic 
adjustment and reduced eustatic sea-level position c. 11.7 ka BP also 
still had a far-reaching impact on the drainage catchments of many 
western European rivers, particularly those feeding towards the 
North Sea. With a land bridge still in place between England and 
continental Europe during the Younger Dryas, the Fleuve Manche in 
the English Channel captured surface drainage from a number of 
major river systems including the Thames, Rhine, and Seine. The 
drainage area of this palaeoriver during the Younger Dryas of c. 
4.97 x 10° km? is estimated to have been much reduced from its 
earlier configuration at the LGM (Fig. 13; Table 2). 


6. Conclusions 


A first-order ice sheet model was used to reconstruct the 
deglaciation of the last Eurasian ice sheet complex (EISC; 23-8 ka 
BP), in a direct continuation of previous work by Patton et al. (2016) 
who modelled the asynchronous growth of the ice complex to its 


Lake volume km? Basin area 10° km? 


= 25.58 
= 13.62 
= 11.96 


— 4.97 
23,464 9.26 
4171 3.61 


Last Glacial Maximum (LGM) extent. Various lines of geophysical 
data were used to guide the broad retreat patterns of the modelled 
ice sheet, including key dating constraints, geomorphological flow 
sets, notable moraine/grounding zone wedge positions, and pat- 
terns of isostatic loading. An evaluation into the robustness of the 
optimum deglaciation experiment against these empirical con- 
straints, as well as a detailed analyses of the results, reveals the 
following insights into the retreat of the EISC: 


e Retreat of the ice complex was asynchronous, reflecting regional 
sensitivities to climate change and various drivers of retreat. 
Even with increases in bulk precipitation, the Celtic ice sheet 
diminished rapidly in response to even modest increases in 
mean annual air temperature, receding to the Scottish High- 
lands by 16 ka BP. However, retreat of the marine-based Barents 
Sea ice sheet was not markedly affected by surface melt, and was 
instead driven by a combination of reduced mass input (pre- 
cipitation) and increased rates of calving. 

e Independent glacial isostatic modelling shows that the pre- 
dicted pattern of ice loading matches reasonably well with Late 
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Glacial relative sea level records from across the Eurasian Arctic, 
particularly over Franz Josef Land, Novaya Zemlya and northern 
Fennoscandia. Model-empirical misfits are best resolved 
through regional optimisations of the Earth model. There is 
support for a low-viscosity zone beneath some regions of the 
Barents Sea, as proposed in earlier work, but spatial differences 
in the optimum Earth model may also reflect shortcomings in 
the ice sheet reconstruction. 

Periodic discharges of icebergs from different regions through 
the Late Weichselian resonate with numerous observations 
from the ice-rafted debris record. These include pronounced 
fluxes at c. 29 ka BP along the western Celtic margin associated 
with early and rapid expansion of ice in this sector, but also the 
collapse of the marine-based Barents Sea ice sheet which 
peaked at 15.2 ka BP with a mass loss rate equivalent to 669 Gt 
a`! The total EISC contribution to Meltwater Pulse 1A, a short- 
lived period of dramatic global eustatic sea-level rise coeval 
with the Bølling warming that exceeded 40 mm per year 
(14.65—14.31 ka BP; Deschamps et al., 2012), was c. 2.5 m of sea- 
level equivalent ice volume. 

Ice caps reacted differently across the domain during the 
Younger Dryas. Increased accumulation of ice across the central 
sectors of the Fennoscandian ice sheet eventually forced a purge 
of ice through the Gulf of Bothnia at 12.3 ka BP. Reduced bulk 
precipitation across the Arctic limits any expansion of the 
remaining ice caps across Svalbard and Russian territories, while 
a cold-based Scottish ice cap grew rapidly from ice-free condi- 
tions within 1 ka. 

The chronology of retreat of Fennoscandian ice through the 
Baltic and Bothnian seas is a major shortcoming of the model 
reconstruction, despite an otherwise sound match to the pattern 
of ice flow history recorded by numerous flowsets. Coupling of 
the ice model to climate modelling output, for example, via 
PMIP4 (Kageyama et al., 2016), may potentially resolve such 
discrepancies in the deglacial mass balance regime, and is an 
area worthy of future inspection. 

The large Baltic and White Sea proglacial lakes expanded rapidly 
after retreat of the FIS during the Bolling/Allerod transition, and 
persisted through the Younger Dryas stadial capturing surface 
runoff from a combined coverage of 12.9 x 10° km? of the 
Eurasian continent at this time. 

Hydrological routing beneath the LGM ice complex reveals 6143 
potential subglacial lake locations across a variety of glaciolog- 
ical settings, with the distribution heavily skewed to surface 
areas < 10 km. The largest lake lies in the Gulf of Bothnia, 
beneath a zone of prominent ice streaming from the Fenno- 
scandian ice sheet, and in combination with a proliferation of 
smaller lakes thus likely formed a major control on its dynamic 
behaviour. 

Catchment areas for the major Eurasian rivers were heavily 
affected by the presence of the EISC as well as eustatic sea-level 
lowering, with the Fleuve Manche river mega-catchment 
dominating the LGM European geography. Its catchment area 
exceeded 2.5 x 10° km? during the Last Glacial Maximum, 
including the drainage of meltwater runoff from one quarter of 
the ice complex. The discharge for such a river system, sensitive 
to periodic meltwater injections triggered by ice sheet fluctua- 
tions, was a probable driver for North Atlantic circulation per- 
turbations for much of the glaciation. 
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